RKKY Interaction in Disordered Graphene 



Hyunyong Lee,^ Junghoon Kim,^ E. R. Mucciolo,^ Georges Bouzerar,^'^ and S. Kettemann^' ^'0 

^Division of Advanced Materials Science, Pohang University of 
Science and Technology (P OS TECH), Pohang 790-784, South Korea 
^Department of Physics, Sungkyunkwan University, Suwon 440-746, Korea 
^Department of Physics, University of Central Florida, Orlando, Florida 32816, USA 
^Institut Neel, CNRS, department MCBT, 25 avenue des Martyrs, BP 166, 38042 Grenoble Cedex 09, France 
^School of Engineering and Science, Jacobs University Bremen, Bremen 28759, Germany 

(Dated: October 31, 2011) 

We investigate the effects of nonmagnetic disorder on the Ruderman-Kittel-Kasuya-Yoshida 
(RKKY) interaction in graphene by studying numerically the Anderson model with on-site and 
hopping disorder on a honeycomb lattice at half filling. We evaluate the strength of the interac- 
tion as a function of the distance R between two magnetic ions, as well as their lattice positions 
and orientations. In the clean limit, we find that the strength of the interaction decays as 1/i?^, 
with its sign and oscillation amplitude showing strong anisotropy. With increasing on-site disorder, 
the mean amplitude decreases exponentially at distances exceeding the elastic mean free path. At 
smaller distances, however, the oscillation amplitude increases strongly and its sign changes on the 
same sublattice for all directions but the armchair direction. For random hopping disorder, no sign 
change is observed. No significant changes to the geometrical average values of the RKKY interac- 
tion are found at small distances, while exponential suppression is observed at distances exceeding 
the localization length. 



An unconventional behavior of the Ruderman-Kittel- 
Kasuya-Yoshida (RKKY) interaction between magnetic 
impurities in undoped graphene was recently reported^ — 
Rather than the conventional decay expected for 

two-dimensional systems, where R is the distance be- 
tween the two magnetic moments, the RKKY interac- 
tion is found to fall off as 1 / R^ at the Dirac (neutrality) 
point. Furthermore, it was found that, due to particle- 
hole symmetry, only ferromagnetic (antiferromagnetic) 
interactions are allowed when two impurities are located 
on the same (different) sublattice.^ 

In a recent experiment, the authors of Ref. measured 
the Kondo effect on graphene samples with a large num- 
ber of vacancies, confirming that these defects induce lo- 
cal magnetic moments.^' - Thus, upon increasing the con- 
trol over the location of such defects, one might be able 
to also measure the RKKY interaction as a function of 
distance and location of local moments. Indeed, a direct 
detection of the RKKY interaction is feasible with the 
recent development of a technique to measure the mag- 
netization curves of individual atoms using spin-polarized 
scanning tunnehng spectroscopy^^*^ With this technique, 
the orientation and distance dependence of the exchange 
interactions can be observed precisely. 

The influence of disorder on the RKKY interaction 
in conventional metals has been thoroughly studied^^ — 
These studies found that the main effect of weak disorder 
is to randomize the electron phase, resulting in an ex- 
ponential decrease of the ensemble- averaged interaction 
amplitude with distance. However, the average does not 
properly characterize the typical interaction strength, as 
any particular disorder configuration has long-range cor- 
relations. Indeed, the typical value, identified as the geo- 
metrical average (J|kky = e^^^/^^ ^"['^^i^ky]'^--^) is found 
to have the same power-law behavior with distance as the 



amplitude of the interaction in the clean limit. Conse- 
quently, at least for conventional metals, weak disorder is 
not likely to cause any critical change in physical proper- 
ties which derive from the RKKY interaction. It is only 
when the system approaches the localized regime that 
the geometrical average is exponentially suppressed^i^ 

In light of these facts, our study focuses on two main 
questions. The first is how a pair of magnetic impurities 
in disordered graphene will interact in general. We con- 
sider impurities located along any lattice orientation, and 
not only along the zigzag and armchair lines. The second 
is how this interaction changes with increasing disorder 
strength. 

Let us begin by considering a general expression for 
the RKKY exchange coupling constant in terms of the 
unperturbed (disorder-free) electronic Green's function 
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Here, J is the local coupling constant between the local- 
ized magnetic impurities and the itinerant electrons, S 
is the magnitude of the impurity spin, i (j) is the site 
index of a magnetic impurity located at position (vj)^ 
f{uj) = [e*^^~^^/^ + 1]~^ is the Fermi-Dirac distribution 
function, and F^^ = tlj^{ri)iljn{rj)tlj^{rj)tljm{ri), with 
ipniXi) denoting the eigenfunction corresponding to the 
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eigenenergy En of the unperturbed electronic Hamilto- 
nian (i.e., in the absence of magnetic disorder). The lat- 
tice constant a and h are set to unity in all calculations. 

Using a zero-temperature approximation (T = 0) and 
changing to an integral form, Eq. ^ can be recast as 
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where F{E,E') = Ke[pji{E)pij{E')\, ji is the Fermi en- 
ergy, and pij{E) = {i\S{E — H)\j)^ which can be calcu- 
lated numerically using the kernel polynomial method 
(KPM).^"^ In the KPM, the matrix elements pij{E) are 
expressed as sums over order-M Chebyshev polynomials 
on the energy E with coefficients obtained through an 
efficient recursion relation involving matrix elements of 
the system Hamiltonian. As our unperturbed electronic 
Hamiltonian with on-site disorder, we employ the single- 
band Anderson tight-binding model on a honeycomb lat- 
tice, 
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where t (~ 2.67 eV for graphene) is the hopping energy, 
Ci {cl) annihilates (creates) an electron at site i, Wi is 
the on-site random disorder energy distributed uniformly 
between [—W/2^ and (i, j) denote nearest-neighbor 

sites. Periodic boundary conditions are used for all cal- 
culations. For clean systems {W = 0), the Chebyshev 
polynomials are calculated up to M = 3 x 10^ on a lat- 
tice with 5 X 10^ sites. 

The RKKY interaction coupling constant between two 
magnetic impurities is calculated using Eq. ^ , of which 
the results for the clean limit are shown in Fig. [TJ In 
order to better visualize the behavior of the amplitude 
in the contour plots, we have multiplied Jrkky by i?^, 
resulting in a smoother (1/i^) decay. The interactions 
along the zigzag and armchair directions are shown sep- 
arately by line plots in Fig. [H These results are in excel- 
lent agreement with previous studiesi — The authors of 
Ref. [3I used a lattice Green's function method to obtain 
an RKKY interaction of the form 
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where all the coefficients are set to unity, = 
(±27r/3A/3, 27r/3) are the Dirac points in the Bloch mo- 
mentum space, R = — r^, and Or is defined in the inset 
of Fig. d^. For a direct comparison, plots of Eqs. (j5j) and 
([6]) are shown in Fig. [1] along with the results calcu- 
lated from Eq. ([3]). As expected from the particle- hole 
symmetry of the spectrum, the magnetic impurity on the 
origin has ferromagnetic correlations with the impurities 
on the same sublattice (Fig.[T^), while antiferromagnetic 
correlations develop for impurities on different sublattices 
(Fig. lib). 
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FIG. 1. (Color online) Plots of the RKKY interaction 
strengths between a magnetic impurity at the origin and an- 
other at: (a) a site from the same sublattice (A A) and (b) 
a site from a different sublattice (AB). In the contour plots, 
the amplitudes are multiplied by the square of the distance 
to facilitate visualization. The lattice constant is set to unity. 
The numerical data is for clean graphene (W = 0). Calcula- 
tions using the kernel polynomial method and lattice Green's 
function method are represented as solid blue and dashed red 
lines, respectively. 



In order to investigate the effect of on-site nonmag- 
netic disorder in graphene, we consider 1.6 x 10^ differ- 
ent disorder configurations for each value of W and then 
evaluate the matrix elements pij through the KPM with 
M = 5 X 10^ on a lattice with 1.8 x 10^ sites. 

For weak (strong) disorder strength, the system is in 
the diffusive (localized) regime, where the actual value of 
W for which this crossover occurs depends on the lattice 
size and has been determined by evaluating the local- 
ization length (see Fig. 2]). The average amplitude of 
the RKKY interaction in the diffusive regime is shown 
in Fig. [21 Similar to conventional metals, the interaction 
decays exponentially with increasing disorder strength as 
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where /g is the mean free path and Jrkky inter- 
action amplitude in the clean limit. It is worth noticing 
that the sign of the interaction oscillates when the impu- 
rities are located along the zigzag- AA direction. 
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FIG. 2. (Color online) Plots of the RKKY interaction 
strength along the (a) zigzag and (b) armchair directions in 
the diffusive regime, as averaged over 1.6 x 10^ different dis- 
order configurations. A lattice with 1.8 x 10^ sites and a 



polynomial degree cutoff of M 
numerical calculations. 



5 X 10 are used in these 



To better characterize the amplitude of the interac- 
tion, we have also calculated the geometrical average 
(^rkky) both diffusive and localized regimes (Fig. [3]). 
In Fig. [3^, one can see that the geometrical average for 
a weakly disordered system remains long ranged and 
has a decaying behavior similar to the clean system. 
As mentioned earlier, studies of conventional metal^ — 
have shown that the geometrical average (i.e., the typical 
value) of the RKKY interaction in weakly disordered sys- 
tems has a power law dependence with the same exponent 
of the clean limit, namely, J^xky ^ 1/^*^ (^-g-, = 2 
in a two-dimensional electron gas). Due to the uncon- 
ventional distance dependence (Eqs. (j5j) and (j6j)) caused 
by the pseudogap at the Dirac point of clean graphene, 
one may expect two possibilities. If the pseudogap is not 
filled by disorder, the geometrical average value of the 
amplitude is expected to have the same exponent of the 
clean system, namely, a = 3. However, if it is filled, then 
the geometrical average value should approach the con- 
ventional power law of a two-dimensional electron gas, 
namely, a = 2. Our calculations show that the former 
is the correct answer. This is in accordance with the 
fact that short-range disorder preserves the pseudogap in 
graphene Therefore, the presence of weak short-range 
disorder in undoped graphene is not anticipated to in- 
duce any major change in physical properties related to 
the RKKY interaction. 

The situation is drastically different in the localized 
regime, where the geometrical average values is exponen- 
tially suppressed with distance, as shown in Fig.[3]3. This 
behavior is captured by the following relation^ 
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where ^ is the localization length. Fig. 2] presents the 
mean free path and the localization length obtained by 
fitting the relations Eq. [3 and Eq. [8] to the numerical 
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FIG. 3. (Color online) Plots of the geometrical average over 
1.6 X 10^ different disorder configurations of the RKKY inter- 
actions for (a) weak and (b) strong disorder. The same lattice 
size and polynomial cutoff of Fig. [2] are used. 



data. For W = t, the localization length is about 10^, 
which is close to the longest linear distance possible in 
our calculations, namely i^max = 60\/3. Therefore, the 
system crosses over from the diffusive to the localized 
regime around W ^ t. For uncorrelated, short-range dis- 
order, which allows for intervalley scattering, the local- 
ization length is given by ^ = /g exp(7rcr/Go)r^i^ where 

^ ~ ^ ('uJa)^+vk4 ' denotes the Fermi velocity, A is 
the energy cutoff, and Go = e'^/h is the conductance 
quantum. 

It is well known that the mean free path is inversely 
proportional to the square of disorder strength (/g ~ 
1/W'^). Therefore, one expects the localization length to 



where 



obey the relation ^ ^ (ci/W^)exp , ^^^^^^^4 
ci is a fitting constant. Indeed, these relations fit reason- 
ably well the numerical data, as shown in Fig. |4j 

To find out the effect of disorder with no sublattice 
symmetry breaking, we added randomness to the hop- 
ping integral and eliminated on-site disorder {wi = 0), 
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where tij = t -\- Atij^ with Atij being distributed uni- 
formly between [—W/2^ W/2]. We perform the same cal- 
culations of the on-site disorder case, but now with a lat- 
tice of 2 X 10^ sites and a Chebyshev polynomial cutoff 
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FIG. 4. (Color online) Plots of the mean free path le and 
the localization length ^ as functions of the disorder strength 
W (in units of t). The blue dashed line represents a fitting 
to the relation le = ci/W^ with a = 70, whereas the red 
dashed line represents the resulting localization length, with 
t'pA = 10 (see text). 
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FIG. 5. Plots of the RKKY interactions along the zigzag- 
AA directions with diagonal disorder (a) W — (h) W = 3t 
and off-diagonal (random hopping) (c) = t, (d) W — 3t 
for 5 X 10^ realizations. A lattice with 2 x 10^ sites and 
a polynomial degree cutoff of M = 10^ are used in these 
numerical calculations. The black dashed line is the averaged 
interaction. 



M = 10^. For comparison, we plot the results together 
with those for the on-site disorder calculations in Fig. [5l 
A total of 5 X 10^ configurations of disorder are used, with 
the thick dashed line indicating the average value. While 
the on-site disorder generates random fluctuations in the 
sign and amplitude of the RKKY interaction (Figs. [5^1, b), 
the hopping disorders affect only the amplitude, even in 
the presence of very strong randomness (Figs. [5]3,d). 



In conclusion, we have confirmed that the RKKY in- 
teractions in clean graphene has a strong anisotropy of its 
sign and oscillation amplitude, and it decays as for 
all directions. Increasing the amount of nonmagnetic, 
on-site disorder causes the averaged amplitude of the 
RKKY interaction to decrease exponentially at distances 
exceeding the elastic mean free path, similarly to what is 
obtained for conventional metals. At smaller distances, 
however, the fluctuations of the amplitude are found to 
increase strongly, with sign oscillations even for a pair 
of magnetic impurities located on the same sublattice, 
for all directions except the armchair direction. When 
the randomness is instead applied to the hopping (off- 
diagonal disorder), the sign oscillations disappear. This 
shows that these sign changes at weak disorder poten- 
tial are caused by the breaking of the sublattice symme- 
try, since off-diagonal disorder preserves this symmetry. 
Our calculations also confirm that the geometrical aver- 
age of the RKKY interaction in disordered graphene has 
the same power law decay at short distances as in the 
clean case. However, it is exponentially suppressed at 
distances exceeding the localization length. We plan, to 
extend these studies by considering the effects of long- 
range disorder and resonant impurities. 
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